Revisit Phase Separation of the Two-Dimensional t — J model by the Power-Lanczos 

Method 
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The power-Lanczos (PL) method is one kind of Green's function Monte-Carlo simulation, which 
is improved by Lanczos iterations. The ground state energies of strongly-correlated models can be 
evaluated by this method quite accurately. In this report, the boundary of phase separation (PS) 
of the two-dimensional t — ,J model is investigated by the power-Lanczos method and Maxwell con- 
struction. The energies are compared with the results evaluated by other methods. Our conclusion 
is that there is no phase separation for ,J/t < 0.4. 



PACS numbers; 74.20.-z, 7L27.-|-a, 74.25.Dw 
I. INTRODUCTION 

It is believed that several key features of the high- 
temperature superconductors (HTSC) can be described 
by the two-dimensional (2D) < — J model on square lat- 
tices. The Hamiltonian is: 

H =-t {C+Cja + h.C.) + J ^ (S, • Sj - jTiiUj), 



<i,j><T 



where < i,j > is the nearest- neighbor pairs and 
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One of the important questions on this model is whether 
the density distribution of charge carriers in the Cu02 
planes is uniform or not. There are several types of 
non-uniform phase such as macroscopic phase separa- 
tion which is separated into a hole- free antiferromagnetir: 
(AF) region and a hole-rich regionEl, and the stripeB tl 
phase which holes form domain walls separating AF re- 
gions. The ground state is determined by the competition 
of the kinetic {t) and exchange (J) terms. The t term fa- 
vors the uniform phase to minimize the kinetic energy, 
while J term tend to attract the electrons together to 
lower the magnetic energy. 

PS state and superconductivity are closely related. 
Some groups argued that the driving mechanism of siu. 
perconductivity is the same as that of phase separationD 
or superconductivity comes from the frustrated phase 
separationQ. On the other hand, if PS occurs in the 
physical regime, the system will become insulating, and 
of course, no superconductivity. Hence it is extremely 
important to determine the phase separation boundary 
of the 2D t — J model to resolve these issues. 

Experimentally, phase separation of the supercon- 
ducting La2Cu04+5 compound are jobserved by neu- 
tron powder [diffractionO, ^'^^La NMRaM and magnetic 
susceptibilityO. La2Cu04+5 phase separates for 0.01 < 



5 < 0.06 below Tp^ « 300K into the nearly stoichiomet- 
ric antiferromagnetic La2Cu04+5j with 5i less than 0.02 
and Neel temperature Tn « 250K, and a metallic su- 
perconducting oxygen-rich phase La2Cu04+52 with 62 ~ 
0.06 and superconducting transition temperature Tc ~ 
34_ftr. In a recent transmission electron microscopy ex- 
periment, an incommensurate modulation was observed 
directly in a phase separated Cu-rich La2Cu04. 003^3. J. 
H. Cho et al£3 measured the magnetic susceptibility of 
the Sr doped compound La2-a;Sr2:Cu04+5. For x < 0.03, 
they also found macroscopic phase separation into super- 
conducting La2_2;Sra;Cu04+5' {S' ~ 0.08) and nonsuper- 
conducting La2-a;Sra;Cu04+6/' {6" « 0.00) phases. For 
X < 0.02, the later phase is antiferromagnetically or- 
dered. The long range order disappears and the phase 
becomes spin-glass-like for 0.02 < x < 0.028. The doped 
holes in the AF (or spin-glass) phase La2_a;Srj;Cu04+5" 
(0 < a: < 0.03, S" fa 0.00) condense into walls separat- 
ing microscopic undoped domains. The phase becomes 
inhomogeneous electronically and magnetically. 

On the other hand, the well known La2-a;Sra;Cu04 su- 
perconductor does not phase separate macroscopically 
because the Sr"*"^ ions are immobile and can not phase 
separate alcpg with the doped holes. In 1989, A. Wei- 
dinger et alS^ observed coexistence of magnetic ordering 
and superconductivity by muon spin rotation (/it+SR) ex- 
periments performed on this material with < a; < 0.15. 
Later D. R. Harshman et a/.li3 interpreted it as the 
phase separation into an antiferromagnetic ordered phase 
and a superconducting phase as in La2Cu04+5 discussed 
above. J. H. Cho et alcB measured the spin dynam- 
ics of this material with 0.02 < x < 0.08 by the ^^^La. 
nuclear quadrupole resonance (NQR) spin-lattice relax- 
ation rates (NSLR) experiments in 1992. They concluded 
that the doped holes were inhomogeneously distributed 
mesoscopically and segregated into walls separating the 
hole-poor antiferromagnetic domains with size « a/^/x « 
10 — lOOA. Later MnS^ ^^^d NQR experiments also sup- 
port this conclusiorOllj. Recent NQR experiment shows 
that in single crystal La2CuO4.02 there are three different 
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regions with different oxygen concentrationglS. 

The measurements of the superconducting dia- 
magnetic moments for underdoped and overdoped 
La2-a:Sr2:Cu04 single crystals show that the underdoped 
sample has only one transition corresponding to Hc2 , and 
the overdoped one has two transitions with the higher one 
at Hc2- The authors proposed that in the two-step tran- 
sition comes from the phase separation in the overdoped 
regimes. 

Several groups have studied the issue of PS by analyt- 
ical or numerical methods. In the low electron density 
region, the theoretical results studied by different meth- 
ods are consistent. Hellberg et al. determined very ac- 
curately that the critical J/t for phase-separation at low 
electron density limit is J/t = 3.4367EI. 

But the theoretical results in the low-dopine_and small 
J/t region are still conflicting. Emery et alh used the 
exact diagonalization (ED) to study the 4x4 cluster. 
The Maxwell construction leads to the conclusion that 
PS ocOjixs. for all values of J/t. Hellberg and Manousakis 
(HM)t3a investigated this problem by the Green Func- 
tion Monte Carlo (GFMC) method and Maxwell con- 
struction for larger clusters and reached the similar con- 
clusion that the t — J model phase separates for all values 
of J/t in the low doping regime. The U{1) slave-boson 
functional integral results reach the conclusion that for 
<j— //i < 0.2 PS occurs for hole density between and 
O.M 

However, the contradictory results are also clairafid by 
several groups. Quantum Monte Carlo (QMC)E3'l3 and 
EEEII results on the Hubbard model, which should be 
consistent with the t—J model for small J/ 1, didn't give 
PS signal. Putikka et al. studied this problem using 
the high-temperature series expansion and found that at 
zero temperature, PjS,|dpes not occur for rJ/i < 1.2 at 
any electron densitynSEa. Prelovsek et aZ.tl used ED to 
study the two-point and four-point density correlations 
on clusters of size 18 and 20 sites. They found that for 
J/t > 1.5 the holes form domain walls along (1,0) or (0,1) 
direction, and phase separate into a hole-rich and a hole- 
free phase for even larger J/t > 2.5. Poilblanc calculated 
the energy of 2 and 4 holes by ED on several clusters up 
to 26 sites. The phase diagram includes a d-wave hole- 
pairing state for J/t > 0.2, a liquid of four-hole droplets 
(quartets) for larger J/rt-|> 0.5, and at even larger J/t, an 
instability towards PSE3. Yokoyama et al. investigated 
the phase diagram by the variational Monte Carlo (VMC) 
methodO. The critical J/t for phase separation at the 
high density limit they found is 1.5, which is consistent 
with Putikka et al. 

The theoretical results of different groups discussed 
above are consistent at the large J/t and low electron 
density region. But in the physical regime of high Tc su- 
perconductors, 0.3 < J/t < 0.5 and high electron density 
0.75 < Ue < 0.95, they are conflicting with eadi. other. 
We have used the power-Lanczos (PL) methodHEj to ob- 
tain the best estimate of the ground state energy in this 
physical regime for the largest cluster (122 sites) that 



have been studied so far. The same-lattices and bound- 
ary conditions as those used by HME3 are studied also in 
this report. It can be shown that our data calculated by 
the PL method, which is a rigorous upper bound of the 
ground state energies, are still lower than extrapolated 
GFMC data reported by HM in the physically interested 
regime. Based on the variational argument we show_that 
there is no phase separation in this physical regimeo In 
this report, we will focus on the competition of the uni- 
form state and PS state. The stripe phase is not taken 
into account. 



II. NUMERICAL METHOD 

The numerical method used here is the "power-Lanczos 
method" , which is a revised vpsion of Green Function 
Monte Carlo (GFMC) methods, and improved by the 
Lanczos iterations such that the ground state properties 
can be calculated more exactly. 

The trial wave functions we used in this report is 
the d-wave-|RVB (resonating valence bond) with AF 
correlationc3, 

I 4-^) = exp{hSM)PdY[{uk + Vkci^clf^^lO) (3) 

with 

"^=2(^"^)'"^^2(^ + ^) 

€k = -2{cosk^ + cosky) - II, Ek = {el + Al) (4) 

where and h are variational parameters controlling 
the magnitude of superconductivity and antiferromag- 
netism. Sm = SI, which is the staggered 

magnetization. operator excludes the states with dou- 
bly occupied sites from the wave function to satisfy the 
constraint of Eq.(^). By tuning the two parameters, the 
variational energy can be minimized and the trial wave 
function corresponding to the best set of parameters will 
be the starting point of the power-Lanczos projection. 

The "power method" is the simplified Green function 
Monte Carlo (GFMC) method, which projects the trial 
wave function determined variationally toward the true 
ground state wave function of the Hamiltonian system- 
atically. 

Suppose the optimal trial wave function is | ^f-^). The 
eigenstates of Hamiltonian are | and Ei are the cor- 
responding eigenenergies. | '5'') denotes the ground state 
wave function. {| ^'')} form a complete set of wave func- 
tions and I '^'^) can be expanded as 

|*^)=^a,|*^) (5) 

i 

Applying the operator {W — H)^ on | 5'"'"), we get 
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i 

= (W^-£;„r^(^^)^a, 1*^) (6) 

i 

where W is some constant. If W is properly chosen such 
that 



W-En \>\ W -E, 



(7) 



for aU i 7^ 0, the weight of excited states will become 
smaller and smaller when the power p is increasing. 
In the p — !■ oo limit, all (-^7^^^)^ become zero and 
I ^1 \I>0). The trial wave function is projected to 

the ground state by applying infinite powers of [W — H). 
For the t — J model, Eq.(^) can be satisfied by setting 
W = 0. If I has the same quantum numbers (S, 
5^, k, etc.) as the ground state | 5'°), and | is not 
orthogonal to | i.e., uq ^ 0, \ = Rp \ 5''^) will 

approach to the ground state as p — > 00. For a physical 
quantity O the expectation value of power p is 



(4r(P) I O I I = 

I HPQHP I *^)/(*^ I H^P I ^'^) 



(8) 



Theoretically we can always get the ground state by 
the power method if we have the trial wave function with 
correct quantum numbers. However, there are two diffi- 
culties with this method. First, the number of interme- 
diate states grow exponentially with powers. Second, the 
statistical error bar grows exponentially with the applied 
power due to the famous 'sign' problem of the Fermionic 
system. We will discuss these two problems below. 

A configuration | a) in real space is defined as 



«) = n <T n 1 0) 



(9) 



where and rj represent the positions of the ith up and 
^th down spins, respectively. \a) satisfies 



H\a) 



J2\a'){a'\H\a) 



(10) 



From Eq.(|), we see that there are three kinds of I a') 
such that the coefficient {a' \ H \ a) 0: flipping a 
nearest-neighbor pair of spin up and down, moving an 
electron from an occupied site to an empty site, and \a) 
itself. The number of such \a') is of order N^, the number 
of electrons. When we apply the next on it, iJ|a') will 
become a summation of about order terms again. To 
measure the expectation value of some physical quantity 
O to the pth power, there are ~ N^p terms to be calcu- 
lated. So it is impossible to calculate physical quantities 
for large power p. To avoid this difficulty, we will not sum 
all the 0{N^P) intermediate states for the chosen config- 
urations. Only the important part of them are taken into 
account by using Monte Carlo method. 



The first step is to generate the sequence of Ntotai con- 
figurations \ tti), I < i < Ntotah by Metropohs algorithm 
which is distributed as the weight P{i) correspondiag 
to the trial wave function optimized by VMC methodEll. 
Our task now is to calculate the matrix elements 



(^^ I HPQHP I 



)o = 



P{^) 



I p 

I HPQHP I ae)a> 



(11) 



where Ci is the amplitude of the configuration | a^) of the 
trial wave function ^I'-^). The operator Q equals to the 
Hamiltonian operator H when we measure the energy. 
The first task is to calculate (i?") which will be useful 
for calculating energy {n = 2p -I- 1), and normalization 
factor (n = 2p). Expanding the matrix element in the 
following form 

(a, |i7"|a/)= J2 {a^\ H \mi){miH \m2) ... 

mi,m2,...,mTi_i 

(m„_2 I H I m„_i)(m„„i \ H \ Uf) (12) 

Define 



Pal3 = P\a}^\P) = -T^Hap 



(13) 



with Za =Y1, pHap and Hap — {a \ H \ f3). Eq.(|T^) can 
be rewritten as 

(fli I iJ" \af) — ^ {Pa )(^mim2^mi) • ■ ■ 

mi .m2 , . ...nin — x 

) (14) 

To calculate this exactly all the possible paths 
{mi, 7712, ■ • ■ , "7„_i} connecting the initial and final con- 
figuration I ai) and | a/) have to be taken into account. 
Since there are about terms, it is impossible to sum 
all the terms. Instead we make use of the random walk 
procedure. 

Starting from the initial configuration | ai), use Paim[ 
defined in Eq. ( p^ ) to be the probability of transition from 
I Ui) to I m'l) state. There are many possible | 777^) 
of three types corresponding to nearest neighbor spin- 
flipping, hole-hopping, and diagonal terms as mentioned 
earlier in this section. For these types the correspond- 
ing matrix elements Ha-m' are —t, — J/2, and — J/2, re- 
spectively. Choose one of these | m[) to be the | nii) 
randomly with the transition probability Paim[- Then 
use the same procedure to generate the next intermedi- 
ate state I 70,2). Repeat this n times, we can arrive the 
final state | a/) of this path. The measurement of (iJ") 
for this path is 



Mi — {Zai ■ Zmi ■ Zm.2 ' ' ' Zm^_-^ )Cj^ / C'^. 



(15) 
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This is only one term in the summation of Eq.(|lJ). Of 
course the error will be large if we use this value to esti- 
mate the summation in Eq.(|lJ). To reduce the error, we 
generate many paths by the procedure above and average 
the measurement 
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Npath r- 



(16) 



Npath is total number of paths of random walk with the 
same initial configuration \ai). {|toJ), |m2, • ■ • , |aj)} is 
the intermediate and final states of the fcth path. Npath 
should be large enough to get good statistics. Typical 
Npath is of order 10^ ~ 10^. Together with Eq.(|ll|), and 
Eq.® we get 



Ntota 



Ntotal 



(17) 



with 



S = 



1 



((Af2) - (Af)2) 



N 



Ntotal 



1 



total 



(18) 



Now we can calculate {H"). With this it is straight- 
forward to calculate the energy of pth power i?^^^ = 

(if2p+l)/(7j2p^. 

Another difficulty is the 'fcrmionic sign problem'. The 
matrix elements of Hamiltionian {a \ H \ P) are not pos- 
itive or negative definite. That is, some of the weight 
Eq.(|l5|) of the Npath paths are positive and others are 
negative. The summation in Eq.([l^), or Eq. (|l7|) will be- 
come a summation of a positive part and a negative part 



{HPORP) = {of + O'^:') ± {aX' + ct'I 



(19) 



where of and O^f ■* are the positive and negative parts 



of the summation in Eq. ([Ij 

For smaller powers, the portion of positive part is much 
larger than the negative part. It is because that a step of 
random walk only flips a pair of nearest neighbor spins 
or moves an electron to the nearest neighbor empty site, 
and the probability of crossing the nodes of the trial wave 
function and resulting a negative matrix element is small. 
For large powers, the probability of crossing the nodes of 
the trial wave function will be accumulated during the 
long path of random walk and the negative part in the 
summation Eq.(p^) and Eq. (p7|) will grow larger. The 
larger the negative portion (absolute value closer to the 
positive portion) the larger the error bars of the physical 
measurements will be. In FigJ^(a) we show the ratio of 
the negative and positive portions varies with power for 
J/t = 0.4 and electron density = 10/16, 10/36, 18/36, 
and 26/36. It is clear that the negative portion is larger 



for higher electron density. For the high density cases 
10/16 and 26/36, the absolute value of negative part is 
almost the same as the positive part for {H"^^). Thus | 
(ff^O) |<| I and the statistical error of (iJ^O) will 

be largely increased. For example, the negative ratio is 
about 79% for 10/36, 87% for 18/36, and 97% for 26/36. 




J'X'l^' O 10/16 



V 18/36 
T 26/36 




FIG. 1. Negative sign ratio which is defined as the ratio of 
the absolute value of negative and positive parts of {H") for 
(a) different electron densities: 10/16 (open circle), 10/36 (full 
circle), 18/36 (open triangle) and 26/36 (full triangle) and (b) 
different J/t values; J/t = (open circle), 0.4 (full circle), 0.8 
(open triangle), 1.2 (full triangle), 1.6 (open square) and 2.0 
(full square). 

If we choose Ntotai and Npath large enough to reduce 
the error bars of (i?^°)± to be as small as 0.1%, the cor- 
responding error bars of of these three electron 
densities will be 1%, 1.5%, and 6.7%. So no matter how 
large Ntotal and Npath are, and how accurate of are, 
the sign problem always prevents us to do calculations 
for large powers, especially in the high electron density 
regime. The largest power we use is usually 6 or 8 for high 
electron densities. The negative ratio is about 80%~90%. 
Beyond this the statistical errors are too large to be reli- 
able. The energies varies with power for different densi- 
ties are plotted in the Fig.|| (open circles) for J/t = 0.4. 
In Fig|^(b) the negative ratio for different J/t value of 
26/36 system is shown. By comparing Figjl](a) and (b) 
we see the negative ratio depends on the electron density 
much more sensitive than J/t. 

The sign problem is the intrinsic limitation of the 
power method. It is almost impossible to project the 
trial wave function to be very close to the ground state 
by applying large powers on it, especially in the physi- 
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cal interested region with high electron density. Other 
methods hke QMC or GFMC also encounter a similar 
problem. To overcome this limitation, we have modi- 
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FIG. 2. Energy vs power for J/t = 0.4 for 10/16, 10/36, 
18/36, and 26/36 (from the top to the bottom, open circles). 
The full circles are the PLl data. 

tied the power method. Here weuwill introduce the more 
powerful power-Lanczos methodEd. 

An inappropriately chosen trial wave function would 
require a lot of computer time to converge to the ground 
state because it is necessary to apply many powers on 
the trial wave function. Hence large Ntotai and Npath 
will be necessary to reduce the statistical error due to 
the random walk approximation and the fermionic nega- 
tive sign problem. Since the negative sign problem pre- 
vents us to do large powers to get the results close to the 
ground state in the physically interested high density re- 
gion, another solution is to make our starting point, the 
trial wave function closer to the ground state. If our trial 
wave function is closer to the ground state, it will not 
be necessary to apply many powers. And the negative 
impact of sign problem can be reduced. Here we use the 
Lanczos algorithm to improve our trial wave function, 
and then apply powers of Hamiltonian on the new trial 
wave function. The method is a combination of power 
method introduced in the previous section and the Lanc- 



zos iteration method usually used in the exact diagonal- 
ization (ED) for the small clusters. Hence it was named 
as 'power-Lanczos method'. Lanczos iteration method 
was first proposed by Heeb ancU-Rice to improve the 
trial wave function systematicalljH. The effectiveness of 
the method is demonstrated on the two-dimensional an- 
tiferromagnetic Heisenherg model and later on the two- 
dimensional t — J modelEEI. We will describe their method 
below. 

Starting from a trial wave function | 0o) (which we 
called it | above) optimized by the variational 

method. Apply the Hamiltonian H on | 0o)- Since | 0o) 
is not an eigenstate of H, it becomes 



with 



I 0o) = ao I 0o) -f- 6i I 01 







(20) 



(21) 



which means | ^o) and | 0i) are orthogonal, oq = (0o I 
H \ 4>o) = Eq, which is the expectation value of energy 
of the trial wave function. And 



bi 



hi 



H 



bo\{H- Eo)H 



(22) 



This formalism is similar to the Lanczos iteration for ex- 
act diagonalization, which the recurrence relation is 

7J I (/.o) = ao I M+bi I (pi) 

H I (j)n) a„ I + bn \ cj)n-l) + &n+l | '/>n+l) (23) 

For small systems, n can be as large as the number of all 
the possible configurations in real space and H can be 
diagonalized exactly in the basis {| 4>q), \ . . . , | 4>n)} 
and the eigenstate corresponding the lowest eigenenergy 
is just the ground state we need. The_largest cluster 
solved by this method exactly is 32 sitesEj. Although for 
larger clusters it is impossible to use this method to solve 
the exact ground state, we can truncate the recurrence 
relation Eq.(p3|) at n which is small enough such that 
we can calculate at large systems, and then diagonalize 
the Hamiltonian by the incomplete basis. The lowest- 
energy eigenstate will not be the ground state because 
the basis is incomplete but it will be much closer to the 
ground state than the original trial wave function even n 
is small. 

Now consider the first order iteration [n — V) which 
we call it 'PLl' in Eq.(i|) 



H\(t>o)=ao \ (t>o)+bi I 4>i) 
H\(l)i)^ai\ (l)i)+bi I 0o) +62 I 02) 



(24) 



where (02 | 4>o) — and (02 J_0i) — 0. Since we truncate 
the recurrence relation Eq.(|2^) at n = 1 here, the last 
term in Eq.(E3) is dropped and the Hamiltonian can be 
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written in the space spanned by the basis {| (/>o), | (f>i)} 
as 



H = 



( Hqo 


Hoi ] 






[ Hw 




\ = 


\ hi ai ) 



(25) 



where Eq — uq. Remember that this is the first order 
approximation because we drop the second order term in 
Eq. (p4[) . Under this approximation, Oi in Eq. (p^) is 

ai = (01 I if I 01) 
1 



'o\{H- Eo)H{H - Eo) 

(00 I (H - Eo)H{H -Eo)\> 
(00 I {H - Eo)H I 0o) 



(26) 



Diagonahze the Hamihonian in Eq. (p5|) , the smaUer 
eigenvalue and the corresponding eigenstate is 



El = i[So + ai - ^{Eo-aif + bfi 



PLl) 



El - ai 



(27) 



We have to know the weight of each configuration in real 
space of the new trial wave function | PLl). If 



PLO) =1 0o) = 



,(0) 



and 



PLl)=^aW|0 



(28) 



(29) 



where af'' is the amplitude of | i). Then 

^-1 , 



= U I PLl) 



(^ I (I PLO) 
, bi 



El ~ ai 
I 0i) 



' Ei-ai 

From Eq.(|2^) we know bi \ (j)i) = {H - Eo) \ PLO). Thus 



(1) (0) 

a.- — a. 



1 



El - ai 
Eo 



{i\{H- Eo) I PLO) 



El - ai 



) + -^Y.^^\H\J)a 
El — ai 



(0) 



Thus the ratio of the weight of | Pil) and | PPG) corre- 
sponding the same configuration | i) is 



,(0) 



Pn 



Pi — ai El — a 



1 a(°) 

(30) 



Heeb and RiceEj proposed to calculate oi and bi by us- 
ing Monte Carlo technique. With ai and 6i, Pi and 



I PPl) can be calculated from Eq.(|2^). However, in this 
method ai and bi have to be calculated very accurately. 
A small error will produce large uncertainty in Pi and 
I PPl). The statistical error is inevitable since they are 
calculated by Monte Carlo method. Here we choose an 
alternative. Since I PPl) in Eq.(|7l) has not been prop- 
erly normalized. All the af'^ can be divided by a constant 
(1 - E^) and Eq.(l^) becomes 



,(0) 



ClY,{^\H\J)- 

3 ' 



(0) 
,(0) 



with 



Ci = 



El — ai— Eo 



(31) 



(32) 



Note that Ci is a constant once the trial wave func- 
tion I PPG) =1 0o) is chosen. Instead of calculating Ci 
by using ai and bi obtained by Monte Carlo method, 
we treat Ci as a variational parameter to optimize the 
energy. The optimized energy is just Pi. For each con- 
figuration the PLO weight of the trial wave function has 
to be calculated first and the PLl weight can be evalu- 
ated by using Eq.(pl]). This procedure largely increases 
the computer time (approximately of order Ne). To sai£e 
computer time Ceperley's algorithm for ratio of weighttll 
is used for calculating a'p /af^ where | i) and | j) are 
connected by P. It is important to check the ratio of 
determinant to prevent the divergence of the ratio. 

and af^ should be calculated directly if their ratio ob- 
tained by Ceperley's method is too large or too small, for 
example, larger than 10^ or less than 10~^. Otherwise 
the results evaluated by using Ceperley's method might 
be in err. 

To estimate the possible range of Ci = l/(-E'i — ai — 
Po), from Eq.(Gfl) it is easy to show 



ai = 



h\{H- Eq)^ I 0o) 
ho\{H- Po)2 I 0o) 



Po 



(33) 



If our trial wave function is not too bad, roughly we can 
estimate that Pi ~ Pq and ai Po. So the range we 
tuning Ci is around —1/Eq. Typical behavior of Pi vs 
Ci is shown in Fig.||. 

From our experience of PLl calculation, the 
=0 energy Pi is approximately equal to the en- 



PLU 



ergy of PLOpo^ 



or PLO 



power— 5 ■ 



Our new trial wave 



function | PLl) is much better than the original trial 
wave function | PLO). The next step is to improve the 
wave function further. There are two ways to achieve this 
purpose: (1) using power method introduced in the pre- 
vious section, applying powers on | PLl) and projecting 
it to the ground state systematically; (2) using the next 
order of Lanczos iteration, i.e., the 'PL2'. 

First let's see the first method: 'PLl-power' calcula- 
tion. We can use the same procedure introduced pre- 
viously to apply HP on the new wave function | PLl). 
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FIG. 3. PLlpoujer=o energy dependence of Ci for 26/36, 
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FIG. 4. PLl negative sign ratio which is defined as the ra- 
tio of the absolute value of negative and positive parts of (H") 
for different electron densities: 10/16 (open circle), 10/36 (full 
circle), 18/36 (open triangle) and 26/36 (full triangle). 



The only difference is that the weight of the configura- 
tions has to be changed from PLO to PLl (Eq.(|3l|)). As 
in the PLO-power calculation, the PLl-power calculation 
is suffered from the fermionic sign problem. In Fig. 4 
the negative sign ratio vs power for | PLl) is shown for 
J/t = 0.4 in 10/16, 10/36, 18/36, and 26/36 lattices. It 
is clear that the behavior of negative sign ratio vs power 
for I PLO) and | PLl) are almost the same, which means 
the maximum power under the limitation of negative sign 
problem is the same for both | PLO) and | PLl). Since 
I PLl) is closer to the ground state than | PLO), the best 
PLl-power energy will be lower than that of PLO. Fig.^ 
shows the PLO and PLl-power energy for J/t = 0.4 in 
10/16, 10/36, 18/36, and 26/36 sites (fuU circles). It is 
clear that PLl always yields lower energies. 

Another way to improve the wave function further is 
to do PL2. Since we have the equations for | PLl), we 
can follow the same procedure which evaluates | PLl) 
from I PLO) to generate | PL2) from | PLl). For a 
configuration | i) the weight corresponding | PL2) is 

af'=a«+C2^(z|i?|j)af) (34) 

3 

which is completely in the same form as Eq. (^T|) . 

Using Eq.(^ll) to replace a'"'^' by the expression of af'^ 
we get 

k 

+ C,Y.{^\H\ j)(af + C\ ^^(j | H \ k)a^^^) 

j k 
(2) (0) 

^ = l + iC\ + C,)J2{^\H\k)^ 

"■i k 



+ C,C,J2i^\H\j){j\H\k)^ (35) 
jk «i 

The I PL2) will be closer to the ground state than | PLl), 
and in principle we can apply powers of H on | PL2) to 
project it toward the ground state. The PL2-powcr re- 
sults should be better than that of PLl-power. With the 
similar rough estimation as for PLl Ci, Ci and|_C'2 in 
Eq.(^5|) are near the value 1/L^o- Chen and Leeo cal- 
culated PL2 in 16 sites. The PL2 energy is very close 
to the ground state energy. In practice, for larger lat- 
tices, it is quite difficult to find the optimal values of Ci 
and 6*2 by using the VMC method because the behavior 
of variational energy in the two-dimensional parameter 
space {Ci,C2} is much more complicated than the 1- 
parameter case in PLl calculation. Another problem is 
that the PL2 calculation costs much more computer time 
than PLO and PLl. We are searching for ways to carry 
out the PL2 calculation more efficiently In this paper, 
most data are evaluated from the PLO and PLl calcula- 
tions. 

Note that the energy calculated from the power- 
Lanczos method is still a variational estimate and an 
upper bound of the ground state energy. 

Here we demonstrate an example of the energies cal- 
culated by the PL method, and compare that with some 
recent data by other groups. For 50 electrons in an 8 x 8 
lattice and J/t — 0.2, the PL energies are shown in TA- 
BLE |. 

It can be seen from TABLE | that the PLl and PL2 
energies converge much faster than PLO (GFMC) ones. 
The upper bound of the ground state energy can be es- 
timated with less powers such that the fermionic sign 
problem is also reduced significantly. 

Another interesting point is the ratio of spin (Es) and 
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kinetic (Ek) energies. Fig.|| shows the two portions of 
energy of different powers for 50/64, J/t = 0.2. It can 
be seen that Eg is overestimated while E^ is underesti- 
mated in the trial wave function. From our experiences, 
Es is usually overestimated for the sake of minimizing 
the energy of the trial wave function in the physical inter- 
ested region of the parameter space. The choice causes 
misjudgements of other physical quantities. For exam- 
ple, the superconducting long range pairing correlation 
function will be strongly overestimated if we choose the 
energy-optimized RVB wave functionO. From the behav- 
ior of Es and Ek , we see that it is a very important task 
to find a new trial wave function having more reasonable 
Eg/Ek value, as well as having lower trial energy. 
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FIG. 5. (a) Es, Ek and (b) Ea/Ek for different powers in 
50/64, J/t = 0.2. 

In conclusion of this section, we proposed a powerful 
method which is a combination of Lanczos iteration and 
power method. The Lanczos iteration is used to improve 
the starting point of the power method which projects 
the new (PL) trial wave function toward the ground state. 
The ground state properties can be estimated more accu- 
rately under the limitation of the fermionic sign problem. 



III. EXTRAPOLATION OF THE ENERGIES 

From TABLE || we conclude that the ground state en- 
ergy of this system isJower than —0.6855 by the vari- 
ational principle. HMEj gives —0.6825(23) by extrapo- 
lating the GFMC data. They used a different method 



TABLE I. Power-Lanczos energies of J/t — 0.2 for 50 elec- 
trons in a 8 X 8 lattice. 



power 


PLO 


PLl 


PL2 





-0.6443(1) 


-0.6709(1) 


-0.6819(11) 


2 


-0.6573(4) 


-0.6768(2) 


-0.6833(11) 


4 


-0.6657(6) 


-0.6803(6) 


-0.6851(29) 


6 


-0.6735(9) 


-0.6820(8) 




8 


-0.6769(12) 


-0.6855(13) 





to calculate energies of all different powers in the same 
time. The energies of different powers calculated by HM 
are close to our PLO-power data. It can be seen that 
the extrapolated GFMC energy is still higher than the 
variational upper bound determined by PLlpower=s or 
Ph2power=4 from TABLE ^. Our energy vs. power curves 
for PLO and PLl are fitted by an exponential function: 

E = Eo + ae-''P (36) 

The results are shown in Fig.0. For PLO data, the fitting 
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FIG. 6. Power-Lanczos energies of 50 electrons in 8 x 8 
lattice, J/t = 0.2. The PLO (PLl, PL2) data are represented 
by open circles (full circles, open triangles). The dotted and 
dashed curves fit the PLO and PLl data, respectively, and 
the horizontal dotted and dashed lines are the extrapolated 
values of power = oo. The horizontal solid line is the ground 
state value estimated by HM. 

parameters are Eq = -0.6881(35), a = 0.0438(33), and 
a = 0.174(26). And Eq = -0.6904(48), a = 0.0193(46), 
and b — 0.158(70) for PLl energies. The most important 
parameter is Eq, the values for PLO and PLl data are 
consistent, which are both much smaller than HM result. 
Since the ground state energy must be lower than the 
lowest power-Lanczos energy, their extrapolated data is 
clearly inaccurate. 

In TABLE ||, we show J/t ^ 0.2 energies for several 
electron densities. The lattices and boundary conditions 
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are chosen exactly the same as those chosen by HeUberg 
et alr3. 

Eextrap are the data extrapolated from GFMG-.and 
Eianc are the quoted values for Lanczos iterationsEj. It 
can be seen that for this J/t, EpL, the PL energies arc 
lower than Eextrap for higher densities, and higher than 
Eextrap for intermediate ones (50/90 and 50/72). 

We emphasize here that in the cases of high density 
and small coupling constant, it is difhcult to get accurate 
enough energies for powers larger than 8 because the sign 
problem enlarges the statistical error bar for such large 
powers. Thus it is unreliable to extrapolate the GFMC 
data to estimate the ground energy in this regime for 
larger clusters, even though the method works well in 
the small clusters like 4 x 4, or at low density and large 
J/t regime because the energies can be calculated very 
accurately for more powers, and the convergence in these 
cases is much faster than the high density, small J/t in 
large clusters. For example, the difference between exact 
energy and the number evaluated by PLl-power=10 for 
10/16 is about only 0.1%. 

TABLE IL Energies of J/t = 0.2 for several densities and 
lattices evaluated by PL method Epl, extrapolated GFMC 
Eextrap, and Lanczos iterations Eiana- 



Ne/Ns 


Epl 


Eextrap 


Elanc 


50/90 


-0.9184(13) 


-0.9246(35) 


-0.9211(54) 


50/72 


-0.8074(11) 


-0.8103(28) 


-0.8098(32) 


50/64 


-0.6855(13) 


-0.6825(23) 


-0.6826(27) 


50/56 


-0.4708(8) 


-0.4681(18) 


-0.4717(23) 


60/64 


-0.3710(6) 


-0.3701(17) 


-0.3693(7) 



IV. DETERMINATION OF PHASE 
SEPARATION: MAXWELL CONSTRUCTION 

There are several criteria to determine the boundary 
of PS. The standard one is to find the density of di- 
vergence of the compressibility (the second derivative of 
the energy-density curve). This method is successful in 
determining the phase boundary of the one-dimensional 
t — J model, which phase separates into one electron-free 
regiop-|and another one containing both electrons and 
holeseil. 

For the 2D t— J model, the PS state contains one hole- 
free Heisenberg AF region and another hole-containing 
one. The finite size effect from the "surface energy" be- 
tween the hole-free Heisenberg region of the PS state is 
significant (of the order 1 / a/ZV^, where Ns is the size of 
the lattice), and this surface energy will push the PS 
state to much larger J/t value. If we use the inverse 
compressibility method in the 2D t — J model, we will 
get wrong results due to this effect. Thus we have to use 
the Maxwell construction method to find the PS bound- 
ary from the energies of the finite-size, uniform systems. 



The energy approaches the ground state as the power 
increases in our PL method. Although in the physical 
regime most of our best data have not yet converged to 
the exact ground state values, the systematic variation of 
energies is enough for us to give a variational estimate of 
the lower bound of the phase separation boundary. We 
will use the lowest PL energies in the following discussion. 
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FIG. 7. The ground state energies of different electron den- 
sities in 50 site lattice for J/t = 2.5. The energies are sub- 
tracted by the perfect PS energy enrie. The straight line 
shows the ground state energy of the PS state by Maxwell 
construction. The onset of PS state is at rie = 0.320. 

Fig.|^ shows the energy-density curve of J/t — 2.5 in 50 
sites. The energies are subtracted by enne, where en is 
the Heisenberg energy and rig the electron density, enne 
is energy of the "perfect phase separation state" , which 
separates into a hole-free Heisenberg antiferromagnet and 
an electron-free vacuum phases. This is the ground state 
of the infinite J/t limit if the quantum fluctuation near 
the boundary of the two regions is not taken into account. 
Thus it is also a good reference state for large J/t cases. 

Fig.|^ shows the PS phase boundary determined by 
Maxwell construction. It is clear that the curve become 
convex for the larger densities and the phase separates 
into a hole-free Heisenberg antiferromagnet and a region 
with density rig = 0.320. The energy of the PS state is 
the linear combination of the two phases, which is repre- 
sented by the jSolid line. The results are consistent with 
Hellberg et aZ.E3, whose critical Ue = 0.296. 

For smaller J/t, it is more difficult to define the tan- 
gent line to find the critical Ue as we did in FigjTj The 
J/t = 0.6 data are shown in Fig.||(a). It can be seen that 
the energy-density curve is almost a straight line for this 
coupling strength. The curve is fitted by a linear func- 
tion, which is shown by the dotted line in Fig.||(a). And 
the deviations of the data points from the line are shown 
in Fig.^(b). We see that the deviations are of the order 
of magnitude of the statistical error bars. Since the sub- 
traction of a linear function won't change the convexity 
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50 sites, J/t=0.6 (a) 
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-0,001 
-0,002 
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(b) 



FIG. 8. (a) Energy-density curve for J/t — 0.6 in 50 sites, 
it can be fitted very well by a straight line, (b) Deviation of 
the energy from the line in (a), the values are of the same 
magnitude of the error bars. 



or concavity of the curve, the Maxwell construction can 
be applied on the curve in Fig.|^(b) as we did in FigJ^. 
A polynomial is used to fit the E — Eunear points. The 
curve is always concave and there is no PS at this J /t. 
An interesting point is that one can always find another 
fitting function gives PS because the E — Eunear value of 
fie = 0.96 is larger than those of Ue = 0.92 and rig = 1. 
Since the E — Eunear values are all so small, the fitting 
is somewhat arbitrary within the error bars. And the 
energy-density curve is almost linear. It is very likely 
that J/t = 0.6 is just the PS phase boundary of vanish- 
ing hole density for 50 sites. 

Another possibility is the formation of stripe state sug- 
gested by S. R. White et alnB. For the stripe state, the 
holes adjust themselves to form stripes to minimize the 
total energy and does not have the problem of the sur- 
face energy as in the PS state. The adjustment makes 
the energy-density curves linear. But since we use the 
uniform trial wave functions and periodic boundary con- 
dition here, it is difficult to get signals of the stripes state. 
This will be our future work. 

From the two previous examples, we see that it is dif- 
ficult to read out the slope variation from the energy- 
density plots to do the Maxwell construction, especially 
for small J/t values. For such J/t, the possible critical 
Tie of PS (if there is) is close to 1. In this regime, the 
energy-density curve is very close to a straight line and 



the ground state energies are very difficult to calculate 
(or estimate) accurately. It will be quite unreliable to es- 
timate the critical Ue and to judge whether there is PS or 
not by extrapolatiott,and curve fitting. Therefore we will 
follow Emery et alH by examining a more well-defined 
value e{x), the energy per hole for hole density x. 

The energy of the PS state can be separated into two 
parts: 



ExN, = {Ns- N)eH + Nen 



(37) 



where Ng is the total number of sites and N is the number 
of sites in the hole-rich phase. The first term of the right- 
hand side of Eq. (|37| ) is the energy of the hole- free Heisen- 
berg region and the second term is that of the hole-rich 
part- en — —1.169 J denotes the Heisenberg energy per 
sitec3. And et is energy per site in the uniform hole-rich 
phase, which is a function of the hole density x — Nh/N. 
Nh is the number of holes. E can be rearranged into the 
form: 



where 



X iV, = NsCH + Nhe(x) 



e{x) = [-en + eh{x)]/x 



(38) 



(39) 



e{x) can be interpreted as the energy per hole relative 
to the half-filled Heisenberg state. If e(x) of a particular 
J/t has a minimum at a; = Xm and the hole density 
of the total system is smaller than Xm, the system will 
adjust the size of the hole- rich phase N such that N^/N is 
equal to Xm and it minimizes the total energy in Eq. (^8|) . 
Since Ns, ch, and Nh. are all constants, the total energy 
is minimized as e{x) is minimized. Thus Xm is the critical 
density for phase separation at this J/t. 

We calculated e(x) from the energy of the uniform 
states eh{x) by the PL method and found the minimum 
of e(x) on several sizes of lattices for several densities 
and J/t. The largest one is the \/122 x \JV12 cluster. 
It is very difficult to get the converged ground state en- 
ergy in the physical regime due to the sign problem. But 
from the systematic PL procedure, we can see the trend 
of the change of e{x) and estimate the boundary of the 
PS state. The PL-1 power=4 (for 82 and 122 sites) or 
PL-1 power=6 (for 50 and 36 sites) energy is used here 
as the eh{x). It is about 2^4 percent lower than the 
variational energy. We estimate the difference between 
the best PL energy is within one or two percentage of 
the true ground state energy. 

In Fig.^we show the e{x) of J/t — 2.5 in 50 sites cor- 
responding to Fig.|^. The data are fitted by a polynomial 
and the minimum of eix) is at a; = 0.68, which is ex- 
actly agree with the density rig = 0.32 determined from 
Fig.0. Thus the two versions of Maxwell construction 
are consistent. The advantage of the latter one is that 
the critical density for PS state will be determined much 
more accurately if the energy-density curve is almost lin- 
ear. This is the case for the physical interested regime. 
So we will use the e(x) to determine the PS boundary in 
this report. 
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FIG. 9. e{x) vs hole density x of J/t = 2.5 in 50 sites. The 
minimum of the fitted curve is at a; = 0.68. 

V. RESULTS AND DISCUSSION 

Using the Maxwell construction we determined the 
phase boundary of PS state for several J/t. The re- 
sults are shown in Fig.|l^. The results evaluated by using 
Green function Monte Carlo with stochastic reconfigu- 
ration (GFMCSR) by M. Calandra et alM and Areen 
function Monte Carlo with extrapolation by HMc2l are 
also shown in the figure. 

From Fig.|l^ we see that in the physical regime, our 
data are consistent with the GFMCSR results that there 
is a lower bound of the J/t value fioj; PS state. This 
is also the case in our previous stud}c3. In contrast, the 
extrapolated GFMC results show that PS state occurs at 
all values of J/t. In the large J/t and low electron density 
region, all the results are qualitatively consistent. 

There are several possible reasons for the inconsistence 
in the physical region. First, the energies calculated are 
not accurate enough. The energy-density curve is very 
close to a straight line in this region, so a small error in 
energy will cause a completely different result of Maxwell 
construction because the method is related to the sec- 
ond derivative of the curve. Fig.|Tl| shows the PS phase 
boundaries determined by energies evaluated by different 
projection levels. It is clear that as the powers increase 
(that is, the energies closer to the ground state energies), 
the PS phase boundary will be pushed to larger J/t value 
and electron densities. This feature appears in all of our 
calculation on different lattices. It is strongly suggestive 
that the PS phase boundary determined from our vari- 
ational calculation is a lower bound of J/t and electron 
densities of the real boundary. 

In TABLE |l| we show that our PL energies are lower 
than the extrapolated GFMC values. And since PL 
method is a variational one, the true ground state ener- 
gies will be even lower. We use J/t — 0.3 in 122 sites as an 
example to demonstrate the effect of the non-converged 
energies in more detail. Again the energy-density plot 




FIG. 10. Phase separation boundary on the phase dia- 
gram of the two-dimensional t — J model evaluated by the 
power-Lanczos method (open circles), GFMCSR method by 
M. Calandra et al. (open squares), and GFMC with extrap- 
olation by C. S. Hellberg et al. (full circles). 

in Fig.^(a) is very close to a straight line and it will be 
dangerous to fit the curve and draw the tangent line from 
the half-filled point. 

In Fig.|l2|(b) we show the e{x) vs. hole density x curves 
for different level of projection. As the power increases, 
the energy decreases and approaches to the ground state 
energy, and e{x) also decreases. It is interesting that e{x) 
of smaller x drops much faster than those of large x. In 
this case, at the PLOpot«er=o level, the minimal e{x) will 
occur at x larger than 0.0656 (8 holes in 122 sites). While 
at the PLlpo^jgT.— Q level, the curve becomes almost flat. 
For PLlpower=4, the smallest density x = 0.0164 (2 holes) 
becomes the lowest one. Note that the energies have not 
yet converged. From the trend of e(x) with increasing 
power, we expect the minimum of e{x) at x = 0.0164 
will be even deeper. That is, the PS state occurs at 
X < 0.0164, while a much larger value xps = 0.123 is 
reported by HM. 

The changes of energies and e{x) with respect to pow- 
ers are shown in Fig.n3. It is clear that e{x) drops faster 
for at lower hole densities with increasing powers. The 
reason is that the hole density is at the denominator of 
the definition of e{x). From the curves in Fig.|l^, it is 
clear that the energies have not converged yet. And the 
trend of the curves strongly suggests that the e{x) of the 
smaller x cases will be much more reduced than those 
of the larger ones. That is, if the ground state energies 
are calculated more accurately, the PS boundary will be 
pushed toward the higher density or larger J/t. Thus 
the true phase boundary will be at the upper right side 
of our variational result in Fig.p^. 

Fig.|lJ shows the similar plots for J/t = 0.4. The be- 
havior is a little subtle in this value of J/t. The e{x) is 
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FIG. 11. PS phase boundaries determined by different PL 
levels. It can be seen that for larger powers, the boundary is 
pushed to larger J jt and electron density rie. 



flatter than the J/t — 0.3 case. The minimal value is 
also aX X = 0.0164 point, but within error bar with the 
X = 0.0328 (4 holes) point. Again the energies are still 
well above the ground state values, we expect the e[x) 
for X — 0.0164 will be even lower. 

Note that in Fig.|l|(b) and Fi^(b) the Heisenberg 
energy used is the thermodynamic limit ch = — 1.169J. 
And en = — 1.1713J for finite system with 122 sites are 
used for Fig.p^(c) and Figp^(c). The later case will make 
the e{x) values larger. And in the J/t — 0.4 case, the 
PLlpou)er=4 curvc becomes almost flat. This flnite-size 
effect will shift the phase boundary a little but the main 
results will be unchanged. 

Another way to understand the finite-size effect is 
to use finite-size scaling. For example, the case of 
J/t = 0.4 in 98 sites shown in Fig. 20 inCJ (calculated 
by M. Calandra et a/.), the minimum of e{x) is at the 
Xm — 0.041 (4 holes) point. It is possible that PS occurs 
at this hole density. But together with the 50 site data 
in the same figure, the Xm — 0.08, also for 4 holes. The 
x„i becomes a half as the size of lattice is doubled. This 
scaling behavior, is consistent with that reported in our 
previous papeiuj. And the scaling beiiavior to even larger 
lattices is also shown in Fig. 3 of Refo. Note that all the 
energies used in these references are upper bounds of the 
ground state energies. We expect that e{x) for smaller 
hole density x will become even lower, and the phase 
boundary of PS state will be pushed to larger electron 
density and largen. J/t. Thus in contrast with the con- 
clusion by in RefEj, we conclude that PS does not occur 
at J/t = 0.4 for any density. 

Another finite size effect comes from the different 
shell for different electron density. The effect makes 
the energy-density curves jagged and difficult to ana- 
lyze. HM calculated the energies for the same number 
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FIG. 12. (a)Energy vs. hole density x for J/t = 0.3 in 
122 sites. (b)PLOpojDei—o (open circle), PLlpomer=o (full cir- 
cle), and Phlpower^A (opcu triangle) e{x) calculated by using 
CH = -1.169J and {c)eH = -1.1713J. 



of electrons (with the same filled shell) in different sizes 
of lattices and different boundary conditions to avoid this 
effect. But there are still other finite size effects in this 
way due to the different sizes and boundary conditions. 
The only way to eliminate the finite size effects is to 
study this problem in different, and large enough sizes of 
lattices. 

For J/t < 0.4, the minimal e{x) is always at two holes 
for lattices of different sizes (up to 1221 This may be 
resulted from the two-hole bound stateEII but not PS at 
Xm — '^/Ng. If there were PS, Xm would be at the same 
(or nearby) density rather than the same number of holes. 

In conclusion, from the analysis of the energies evalu- 
ated from the PL method, we conclude that J/t — 0.4 
is a lower bound for PS. The phases will not separate 
for J/t smaller than this value. From the trend of the 
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FIG. 13. (a) e(x) of J/t = 0.3 for different powers in 122 
sites. The electron numbers are 120 (open circles), 118 (full 
circles), 116 (open triangles) , and 114 (full triangles). (b)Tlie 
ratios of power energies EpL and variational {PLOpomer=o) 
energy Eo- (c)The differences between e{x) of different powers 
and the variational value. 



energies and e{x) functions, this lower bound may be 
pushed to even larger J/t. The extrapolation of ground 
state energies and energy-density curves may be unreli- 
able and misleading in the physical interesting regime. F. 
Becca et al. give upper and lower bounds for the ground 
state energy of the infinite-C/ 2D Hubbard model, which 
is equivalent to the J/t = case for the t — J model. 
Their analysis on this case ruled out the possibility of 
PS for the electron density less than half-filleal3. By the 
way, from the finite size analysis of the exact diagonaliza- 
tion data for one and two holes in smaller clusters (up to 
32 sites) gives the reaulUJ^hat there is no two-hole bound 
state for J/t < 0.6OEZlca. Thus it seems the attractive 
interaction is not strong enough to be the mechanism to 
cause the PS state for small J/t values. 
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